#define MAX_SIZE 1111

*****************************************************************************
*									    *
*	1D (slow) Fourier Transform					    *
*									    *
*****************************************************************************
      subroutine dft1du( plusminus, n, a, inc)
      double precision a(*)
      double precision w(MAX_SIZE)
      integer n, plusminus, inc
      
      double precision u, v, tpi
      TPI = 2 * 3.14159265358979323846

      if( inc .ne. 1) goto 11
      if ( plusminus .eq. -1) then
	l = (n+2)/2
c$doacross local(k,u,v)
	do k = 1, l
	    u = 0.0
	    v = 0.0
	    do i = 1, n
		u = u + a(i) * cos((k-1) * (i-1) * tpi / float(n))
		v = v - a(i) * sin((k-1) * (i-1) * tpi / float(n))
	    end do
	    w((2*k-2)+1) = u
	    w((2*k-1)+1) = v
	end do
c$doacross
	do i = 1, 2*l
	    a(i) = w(i)
	end do
      elseif (plusminus .eq. 1) then
	Print *,'NOT Implemented'
      else
	print *, 'Error first argument of DFT1DU should -1 or +1'
	stop
 11	print *, 'This Simple version support only INC = 1'
	stop
      endif

      return
      end


      subroutine zft1d( plusminus, n, a, inc)
      implicit double precision ( a-h, o-z)
      double complex a(*)
      double complex w(MAX_SIZE)
      integer n, plusminus
      
      double complex cp
      double precision re, im, tpi

      tpi = 2.*3.14159265358979323846

	do i = 1, n
	    w(i) = 0.0
	    do k = 1, n
		im = mod(((i-1)*(k-1)), n)
		im = dfloat( plusminus*im)*tpi / dfloat(n) 
		re = cos( im)
		im = sin( im)
		cp = dcmplx( re, im)
		w(i) = w(i) + cp*a(inc*(k-1)+1)
	    end do
	end do
	do i = 1, n
	    a(inc*(i-1)+1) = w(i)
	end do

      return
      end



*****************************************************************************
*									    *
*	2D (slow) Fourier Transform					    *
*									    *
*****************************************************************************
      subroutine zft2d( plusminus, n1, n2, w, ldw)
      double complex w(ldw,n2)
      integer n1, n2, lda, ldw, plusminus

c$doacross local(j)
      do j = 1, n2
	call zft1d( plusminus, n1, w(1,j), 1)
      end do

c$doacross local(i)
      do i = 1, n1
	call zft1d( plusminus, n2, w(i,1), ldw)
      end do

      return
      end

      subroutine dft2du( plusminus, n1, n2, w, ldw)
      Double Precision w(ldw,n2)
      integer n1, n2, lda, ldw, plusminus, l, lh
 
      l = ((n1+2)/2)*2
      lh = ldw/2
      if (ldw .lt. l) then
	PRINT *, 'ERROR Leading Dimension is not sufficient',short(l),short(ldw)
      endif
      if( mod(ldw,2) .ne. 0) then
	PRINT *, 'Warning ... this Checking routine assumes an even leading Dimension '
      endif 
      if (plusminus .eq. -1) then
c$doacross local(j)
	do j = 1, n2
	    call dft1du( plusminus, n1, w(1,j), 1)
	end do
c$doacross local(i,j,my)
	do i = 1, l, 2
	    call zft1d( plusminus, n2, w(i,1), lh)
	end do
      else
c$doacross local(i,j,my)
	do i = 1, l, 2
	    call zft1d( plusminus, n2, w(i,1), lh)
	end do
c$doacross local(j)
	do j = 1, n2
	    call dft1du( plusminus, n1, w(1,j), 1)
	end do
      end if
      return
      end
      

*****************************************************************************
*
*	3D (slow) Fourier Transform
*
*****************************************************************************

      subroutine zft3d(plusminus,n1,n2,n3,w,ld1,ld2)
      double complex w(ld1,ld2,n3)
      integer n1,n2,n3,ld1,ld2,plusminus
c
c   transform along X  first ...
c
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call zft1d( plusminus, n1, w(1,j,k), 1)
       end do
      end do
c
c   transform along Y then ...
c
c$doacross local(i,j,k)
      do k = 1,n3
       do i = 1,n1
	    call zft1d( plusminus, n2, w(i,1,k), ld1)
       end do
      end do
c
c   transform along Z finally ...
c
c$doacross local(i,j,k)
      do i = 1, n1
       do j = 1, n2
	    call zft1d( plusminus, n3, w(i,j,1), ld1*ld2)
       end do
      end do

      return
      end

      subroutine dft3du(plusminus,n1,n2,n3,w,ld1,ld2)
      Double Precision w(ld1,ld2,n3)
      integer n1,n2,n3,ld1,ld2,plusminus, l

      if( plusminus .ne. -1) then
	print *, 'Option non implemented YET'
      end if
      if( mod(ld1,2) .ne. 0) then
	PRINT *, 'Warning ... this Checking routine assumes an even leading Dimension '
      endif 
      l = ((n1+2)/2)*2
      if( ld1 .lt. l ) then
	print *, 'Leading dimension not sufficient'
      endif
c
c   trandform along X  first ...
c
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call dft1du( plusminus, n1, w(1,j,k), 1)
       end do
      end do
c
c   trandform along Y then ...
c
c$doacross local(i,k)
      do k = 1,n3
	do i = 1, l, 2
	    call zft1d( plusminus, n2, w(i,1,k), ld1/2)
	end do
      end do
c
c   trandform along Z finally ...
c
c$doacross local(i,j)
      do j = 1,n2
	do i = 1, l, 2
	    call zft1d( plusminus, n3, w(i,j,1), (ld1*ld2)/2)
	end do
      end do

      return
      end

*****************************************************************************
*									    *
*	1D (slow) Fourier Transform					    *
*									    *
*****************************************************************************
      subroutine sft1du( plusminus, n, a, inc)
      real a(*)
      real w(MAX_SIZE)
      integer n, plusminus, inc
      
      real u, v, tpi
      TPI = 2 * 3.14159265358979323846

      if( inc .ne . 1) goto 11
      if ( plusminus .eq. -1) then
	l = (n+2)/2
c$doacross local(k,u,v)
	do k = 1, l
	    u = 0.0
	    v = 0.0
	    do i = 1, n
		u = u + a(i) * cos((k-1) * (i-1) * tpi / float(n))
		v = v - a(i) * sin((k-1) * (i-1) * tpi / float(n))
	    end do
	    w((2*k-2)+1) = u
	    w((2*k-1)+1) = v
	end do
c$doacross
	do i = 1, 2*l
	    a(i) = w(i)
	end do
      elseif (plusminus .eq. 1) then
	print *,'NOT Implemented'
      else
	print *, 'Error first argument of DFT1DU should -1 or +1'
	stop
 11	print *, 'This Simple version support only INC = 1'
	stop
      endif

      return
      end


      subroutine cft1d( plusminus, n, a, inc)
      implicit real ( a-h, o-z)
      complex a(*)
      complex w(MAX_SIZE)
      integer n, plusminus
      
      complex cp
      real re, im, tpi

      tpi = 2.*3.14159265358979323846

	do i = 1, n
	    w(i) = 0.0
	    do k = 1, n
		im = mod(((i-1)*(k-1)), n)
		im = float( plusminus*im)*tpi / float(n) 
		re = cos( im)
		im = sin( im)
		cp = cmplx( re, im)
		w(i) = w(i) + cp*a(inc*(k-1)+1)
	    end do
	end do
	do i = 1, n
	    a(inc*(i-1)+1) = w(i)
	end do

      return
      end



*****************************************************************************
*									    *
*	2D (slow) Fourier Transform					    *
*									    *
*****************************************************************************
      subroutine cft2d( plusminus, n1, n2, w, ldw)
      complex w(ldw,n2)
      integer n1, n2, lda, ldw, plusminus

c$doacross local(j)
      do j = 1, n2
	call cft1d( plusminus, n1, w(1,j), 1)
      end do

c$doacross 
      do i = 1, n1
	call cft1d( plusminus, n2, w(i,1), ldw)
      end do

      return
      end

      subroutine sft2du( plusminus, n1, n2, w, ldw)
      real w(ldw,n2)
      integer n1, n2, lda, ldw, plusminus, l, lh
 
      l = ((n1+2)/2)*2
      lh = ldw/2
      if (ldw .lt. l) then
	PRINT *, 'ERROR Leading Dimension is not sufficient',short(l),short(ldw)
      endif
      if( mod(ldw,2) .ne. 0) then
	PRINT *, 'Warning ... this Checking routine assumes an even leading Dimension '
      endif 
      if (plusminus .eq. -1) then
c$doacross local(j)
	do j = 1, n2
	    call sft1du( plusminus, n1, w(1,j), 1)
	end do

c$doacross local(i,j,my)
	do i = 1, l, 2
	    call cft1d( plusminus, n2, w(i,1), lh)
	end do
      else
c$doacross local(i,j,my)
	do i = 1, l, 2
	    call cft1d( plusminus, n2, w(i,1), lh)
	end do
c$doacross local(j)
	do j = 1, n2
	    call sft1du( plusminus, n1, w(1,j), 1)
	end do
      end if
      return
      end
      

*****************************************************************************
*
*	3D (slow) Fourier Transform
*
*****************************************************************************

      subroutine cft3d(plusminus,n1,n2,n3,w,ld1,ld2)
      complex w(ld1,ld2,n3)
      integer n1,n2,n3,ld1,ld2,plusminus
c
c   transform along X  first ...
c
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call cft1d( plusminus, n1, w(1,j,k), 1)
       end do
      end do
c
c   transform along Y then ...
c
c$doacross local(i,j,k)
      do k = 1,n3
       do i = 1,n1
	    call cft1d( plusminus, n2, w(i,1,k), ld1)
       end do
      end do
c
c   transform along Z finally ...
c
c$doacross local(i,j,k)
      do i = 1, n1
       do j = 1, n2
	    call cft1d( plusminus, n3, w(i,j,1), ld1*ld2)
       end do
      end do

      return
      end

      subroutine sft3du(plusminus,n1,n2,n3,w,ld1,ld2)
      real w(ld1,ld2,n3)
      integer n1,n2,n3,ld1,ld2,plusminus, l

      if( plusminus .ne. -1) then
	print *, 'Option non implemented YET'
      end if
      if( mod(ld1,2) .ne. 0) then
	PRINT *, 'Warning ... this Checking routine assumes an even leading Dimension '
      endif 
      l = ((n1+2)/2)*2
      if( ld1 .lt. l ) then
	print *, 'Leading dimension not sufficient'
      endif
c
c   trandform along X  first ...
c
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call sft1du( plusminus, n1, w(1,j,k), 1)
       end do
      end do
c
c   transform along Y then ...
c
c$doacross local(i,k)
      do k = 1,n3
	do i = 1, l, 2
	    call cft1d( plusminus, n2, w(i,1,k), ld1/2)
	end do
      end do
c
c   transform along Z finally ...
c
c$doacross local(i,j)
      do j = 1,n2
	do i = 1, l, 2
	    call cft1d( plusminus, n3, w(i,j,1), (ld1*ld2)/2)
	end do
      end do

      return
      end

